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Abstract 

The critical behavior of charge-density waves (CDWs) in the pinned phase is 
studied for apphed fields increasing toward the threshold field, using recently 
developed renormalization group techniques and simulations of automaton 
models. Despite the existence of many metastable states in the pinned state 
of the CDW, the renormalization group treatment can be used successfully 
to find the divergences in the polarization and the correlation length, and, to 
first order in an e = 4 — d expansion, the diverging time scale. The automa- 
ton models studied are a charge-density wave model and a "sandpile" model 
with periodic boundary conditions; these models are found to have the same 
critical behavior, associated with diverging avalanche sizes. The numerical 
results for the polarization and the diverging length and time scales in dimen- 
sions d = 2,3 are in agreement with the analytical treatment. These results 
clarify the connections between the behaviour above and below threshold: 
the characteristic correlation lengths on both sides of the transition diverge 
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with different exponents. Tfie scaling of the distribution of avalanches on the 
approach to threshold is found to be different for automaton and continuous- 
variable models. 
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I. INTRODUCTION 



The static and dynamic behaviour of shding charge density waves (CDWs) is perhaps 
the most well studied example of a class of problems involving the transport of an elastic 
medium through a disordered background. The CDW, which behaves like an elastic medium, 
is pinned by impurities distributed randomly throughout the material. As the magnitude of 
an externally applied electric field is varied, a depinning transition is seen, from a stationary 
phase at weak fields, to a moving phase at strong fields where the CDW slides through the 
material. In the vicinity of this depinning transition, the dynamics of the CDW are correlated 
over long distances, with characteristic correlation lengths diverging at the threshold field. 
It has been shownlll that the behaviour near the threshold field can be studied as a critical 
phenomenon associated with a second order phase transition.! 

Extensive numerical simulationsi^ on a number of classical models for CDWs have 
helped in understanding the critical properties in the vicinity of the depinning transition. 
In particular, critical exponents describing the scaling of physical quantities can be deter- 
mined for various spatial dimensions of the CDW. Recently, it has also proved possible to 
obtain analytically the behaviour in the moving phase above threshold, through a renormal- 
ization group treatment.il The results, obtained within an e-expansion for = 4 — e spatial 
dimensions, agree very well with the numerical simulations. 

While the critical behaviour above threshold is fairly well understood, it is more difficult 
to carry out an analytical treatment below threshold. The reason why the analysis is simpler 
in the moving phase is that, for models in which the elastic interactions within the CDW 
are strictly convex, and no dislocations ( "phase slips" ) are allowed, the system approaches a 
unique steady-state configuration at long times (up to time translations) .0 A perturbation 
expansion around some relatively simple configuration is thus more likely to succeed. In 
the stationary phase, on the other hand, the impurities can pin the CDW in any of several 
different states; the particular state that the CDW reaches depends on the past history of 
the system. Any analytical treatment must include this history dependence, instead of being 
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the simple search for a single attractor, as is sufficient above threshold. 

The history dependence of the CDW below threshold is most easily seen in the polariz- 
ability. If the driving force F (produced by the external electric field) is raised monotonically 
towards the depinning threshold Ft-, the response x^{F) to infinitesimal increases in the force 
becomes larger and larger, and diverges at Ft as x\F) ~ {Ft ~ F)~"' , with an exponent 
7. Also, at any value of F, the response of the CDW to increasing the force further is 
dominated by localized regions of activity, or "avalanches" , whose characteristic size can be 
used to define a diverging correlation length ^ ~ {Ft — F)^'^ . If instead of this monotonic 
approach, the force is lowered from just below threshold, the polarizability x^{F) does not 
diverge at Ft (although, depending on the specific model used for the CDW, it may have 
a cusp singularity). It is thus necessary to distinguish between the upwards polarizability, 
X^(-F), and the downwards polarizability x^{F)- 

In this paper, we show how it is possible to obtain analytically the critical behaviour 
for the 'natural' monotonic approach to threshold, through an expansion around the mean 
field solution with the proper history dependence. Our results are in good agreement with 
recent numerical simulations, as well as new data that we present here. We find that the 
correlation length exponent is given by 

V = 2/d (1) 

for a d dimensional system, while the polarizability exponent is 

7 = (2) 

These exponents can in fact be obtained for a highly simplified CDW mo del0 that can be 
solved exactly. However, the dynamic exponent z cannot be obtained from this model: we 
show that z is equal to its value above threshold 

z = 2-e/3 + 0{e^) (3) 

in d = A — e dimensions. Later in this paper, we discuss how the model of Ref. |TT] corresponds 
to a pathological limit, and why, despite this, it obtains u and 7 correctly. 
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The solution that we obtain is similar in form to the solution above threshold. □ However, 
there are complications with two-sided scaling. This is because, as was seen above threshold, 
there is an extra relevant operator at the fixed point of the renormalization group, which 
affects only the statics and not the dynamics. This operator has very different consequences 
in the pinned and moving phases of the CDW, resulting in different correlation length 
exponents on the two sides of the phase transition. The absence of two-sided scaling is unlike 
the behaviour for phase transitions in equilibrium systems, and conventional arguments 
about the properties of scaling functions must be used with caution here. 

In "phase-only" models for CDWs, the CDW is described by a periodic modulation in the 
density of electrons in the material; the phase of this modulation is assumed to be the only 
dynamical variable of importance and is assumed to be a continuous function of position. 
In numerical simulations, the system is discretized to a lattice, with a continuous phase 
variable at each site. Although simulations on such continuous variable models have been 
performedS^, it is much more efficient to simplify the model, with the phase of the CDW 
at any lattice site restricted to being an integer multiple of 2tt away from a (site dependent) 
preferred phase, and time being incremented in discrete steps.i'i The resulting automaton 
models, which differ in details of how these approximations are made, yield far more precise 
results. Within the numerical uncertainties, much of the critical behaviour is the same for 
both classes of models. In particular, the correlation length exponent u and (with larger 
numerical uncertainties) the polarizability exponent 7 appear to be the same. An important 
exception to this is the downward polarizability, x^'- ^is already mentioned, the nature (or 
even the existence) of a cusp singularity in does depend on the details of the model used, 
and is in fact even different in different continuous variable models. However, this difference 
is not seen in the monotonic approach, which is what we shall be interested in here. 

Since the approach to threshold causes the system to become increasingly unstable, with 
larger and larger avalanches and a divergent polarizability, it is perhaps not surprising that 
differences between automaton and continuous dynamics do not affect the critical behaviour. 
An important exception to this rule is the distribution for the number of avalanches. For 
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continuous dynamics, an infinitesimal increase in F triggers avalanches in which the change 
in the phase at any point is bounded above by 27r (corresponding to unity in the automaton 
models). Most of the avalanche area advances by almost 27r. In large avalanches, it is highly 
probable that a subsequent small increase in F will "retrigger" avalanches in the same region, 
as the original unstable point is likely to be nearly unstable at the completion of the first 
avalanche. Avalanches of usually decreasing size will be retriggered until the region is more 
stable. In the automaton models, owing to the discretization of the CDW phases, several 
"retriggered" avalanches are grouped together into a single large avalanche. This difference 
between the two classes of models is one that involves the behaviour at short times] this 
does not alter long-wavelength low-frequency characteristics that determine v and z (nor 7, 
which involves a spatial average over the entire system). However, it may affect avalanche 
distributions, which do depend on how avalanches are grouped. In this paper we present 
numerical data on the distribution of avalanches for the automaton models, and examine 
the differences with results for continuous variable models. 

The automaton models for CDWs are closely related to some of the sandpile models that 
have been proposed and studied recently.0 We find that the avalanche distribution that we 
obtain here for the CDW automaton models agrees very well with the numerical results on 
sandpiles at the critical point .0 The dynamic exponent, z, that relates the characteristic 
duration of an avalanche to its size, is also found to be the same numerically for both classes 
of models: 

z{d = 2) = 1.32 ±0.04 

z{d = 3) = 1.65 ±0.06 (4) 

The approach to threshold for sandpile models has been treated through numerics and scaling 
arguments by Tang and Bak0; our results differ in certain respects from theirs, which we 
discuss in detail later in this paper. 

The connection between sandpiles and CDW models can be exploited to obtain the 
dynamic exponent z for two dimensional systems. Majumdar and DhaiS have shown that 
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this exponent is equal to 5/4 for sandpile models in d = 2; this is in fair agreement with 
the numerical results for CDWs.lll They are also able to relate the two exponents that enter 
the avalanche distribution to each other, thus leaving only one of the two to be determined 
numerically. In this paper we show that it is in fact possible to obtain the same exponent 
identity from completely different scaling arguments on CDWs, if we assume that the total 
rate of avalanche generation (which is dominated by small avalanches) is not singular as 
threshold is approached. (We have verified this assumption numerically for the automaton 
models.) 

The rest of this paper is organized in the following manner. Section II carries out the 
analytical treatment of the behaviour below threshold, and obtains the critical exponents. 
Section III examines the numerical results, which agree quite well with the analytical predic- 
tions. Section IV compares the distribution of avalanches for CDW continuous dynamics, 
CDW automata, and sandpile models, as well as exploring other connections to sandpile 
models. Section V discusses the relationship of this work to earlier results, and its possible 
relevance to other physical systems. 

II. ANALYTICAL RESULTS. 

The model that we use in this section is the Fukuyama- Lee- Rice mo del,§ which assumes 
that the distortions of the CDW are continuous {i.e. that there are no dislocations). The 
dynamics can then be expressed in terms of a phase variable, (j){x]t), which measures the 
distortions with respect to an ideal undistorted CDW. (The position x is a d-dimensional 
vector.) Assuming that the dynamics are strongly overdamped, and are given by a simple 
relaxation of an energy functional H, the equation of motion is 

d(j){x; t)/dt = -dH/d(j) = V^(/)(x; t) + F + h{x)Y{(j) - f3{x)). (5) 

In Eq.(^) there are three terms on the right hand side: a simple elastic force that arises 
from an elastic energy i(V0)^, a uniform force F from the external electric field, and a force 
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h{x)Y {(f) — f3{x)) from the impurity pinning. This impurity force has an exphcit dependence 
on the position x, arising from the quenched randomness in the location of the impurities, 
through the functions h{x) and These correspond physically to the strength of the 

impurity pinning, h{x), and a preferred phase P{x) selected by the impurities. (In order 
to fix the normalization of h, we choose \Y\ to have a maximum value of 1.) Since the 
impurities have only short range correlations, h{x) and f3{x) are taken to be uncorrelated 
from one position to another, with a distribution p{h) for h and a uniform distribution 
over [0, 2-7?) for f3. The pinning force is the derivative, —d^V{(f); x) of an impurity pinning 
potential. Because of the 27r translational invariance of the CDW, the force Y has to be 
27r-periodic in 0. 

Among the quantities of physical interest in this model below threshold are the polar- 
ization, P, defined as the mean displacement from some initial configuration of the phase 0; 
close to threshold, the singular part of P scales as 

Ps{F) = m - (0)initiai]. ~ {Ft - F)-^+\ (6) 

The derivative of this with respect to F is the upwards polarizability, defined earlier as the 
response to an infinitesimal increase in F, which scales as 

X^F) = ^lim^^[P(F + AF) - P{F)]/AF ~ {Ft - F)-\ (7) 

In Ref. the behaviour above threshold for models of this type was analyzed. This 
was done by expanding around the solution of Eq.(^ within a mean field approximation, 
which is obtained by replacing the short-ranged elastic term with an infinite-ranged 
interaction, (0) — 0. The mean field solution is known to depend on the details of the 
pinning potential V{(f)): for instance, the velocity above threshold scales as f ~ (F — Ft)^ 
in the critical regime,^ with an exponent (3mf = 3/2 for smooth pinning potentialsS and 
/9a/f = 1 for pinning potentials which have wedge-shaped linear cusped maxima.i This lack 
of universality arises from the fact that, for smooth pinning potentials, when the phase 
at any site passes through an unstable point in its local effective potential, it spends a long 



time accelerating before jumping forward to a new stable state. This results in a second time 
scale that diverges at threshold (in addition to the natural scale of lT\jv). No such second 
divergent time scale is present for linear cusped potentials, for which the acceleration time 
remains of 0(1) arbitrarily close to threshold. Despite this difference in behaviour for the two 
types of potentials in mean field theory, numerically it is found that the critical behaviour 
is independent of the shape of the pinning potential for < 4, and is in fact even the same 
for automaton models. In Ref. P, it was argued that this is because the dynamics for small 
d are very irregular: when the phase at any site jumps forward, it produces a sharp force on 
all the neighbours to which it is elastically coupled. For large rf, the total elastic force acting 
on any phase is the sum of contributions from a large number of neighbours; however, as d is 
decreased, the elastic force becomes more and more jerky. For sufficiently small d, the phase 
of the CDW at any site is 'kicked' by its neighbours over the maximum of the local effective 
potential it experiences. The details of the pinning potential are no longer important, and 
the second divergent time scale is eliminated, restoring universality. The expansion around 
mean field theory was accordingly carried out for linear cusped potentials, for which the 
second time scale is absent initially, so that the expansion is better behaved. 

When Ft is approached monotonically from below threshold, the critical behaviour is 
controlled by large avalanches, in which the CDW phase in localized regions moves forward 
abruptly in response to a small change in the field. Since the collapse to universality above 
threshold was understood in terms of the jerky dynamics, the same should be true even below 
threshold. This is indeed found to be the case numerically. Thus even below threshold, we 
seek to expand around the mean field solution for linear cusped potentials. 

Before proceeding with the calculations, it is necessary to consider a possibility that might 
make this procedure erroneous. While we have answered the question of history dependence 
by choosing a specific approach to threshold, namely increasing F monotonically, this does 
not completely eliminate problems that arise from the existence of many metastable states. 
An expansion around the mean field solution of the type we develop here corresponds to first 
increasing the force within an infinite ranged model, and then tuning the elastic coupling to 
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a short ranged form {i.e., restoring the term in Eq.(^ which was replaced by (0) — in 
mean field theory). On the other hand, the physically relevant critical behaviour corresponds 
to first tuning the elastic coupling to a short ranged form, and then increasing the force 
within a short ranged model. When there are many metastable solutions, as is the case 
below threshold, the two methods could in principle yield different critical behaviour. 

Since the CDW phase at any site moves only forward when the force is raised mono- 
tonically, such problems from metastability could arise if reducing the range of the elastic 
interactions from a mean field theory were to result in the phases at different sites trying 
to move backwards. While this will of course be true for some of the phases, for which the 
coupling to their neighbours tends to retard their motion, we shall see later in this paper 
that the polarizability exponent 7 is in fact larger for d < 4 than in mean field theory. Thus 
the dominant effect of reducing the range of the elastic interactions is to make the phases 
move forward, making the polarization (and hence the polarizability) more divergent than 
in mean field theory. We therefore expect our procedure of expanding around the mean field 
solution to yield the correct critical behaviour. 

The expansion around mean field theory is carried out following the prescription of 
Sompolinsky and Zippelius;^ the details are given in Ref. |^. The method yields an impurity 
averaged generating functional for the correlation and response functions of the phase 0, of 
the form 

Z = J[d<!?] [ d^exp{- J d'^xdt^{x,t)[dt-V^ + r](^{x,t) + ^{x,t)[F - FMF{{(l){t)))] 

- J d'^xdtidt2 ti)^{x, t2)C[(0(ti)) + <l>(x, ti) - (0(t2)) - ta)]}. (8) 

Here $ and $ are dummy field variables, in terms of which the long wavelength low frequency 
forms of the truncated correlation and response functions of are given by 

d''{(f){xi,ti) . . .0(x„ , t„)) 

trunc 

/de{ 

= {^{Xi, ti) ... ^{Xm, tm)^{Xm+l,tm+l)^{Xm+n, tm+n)) (9) 

where the left hand side represents the generalized response to a perturbing force e{x,t) 
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added to the right hand side of Eq.(^, and the right hand side is the expectation value 
using the measure of Eq.(|^). is hke a coarse grained version of the phase (j){x,t). It 

gives the deviations of 0(a;,t) from its spatially averaged value, (0), and must satisfy the 
consistency condition ($) = 0, which determines F — Fmf- Various other terms, that do 
not affect the critical behaviour, have been suppressed in Eq.(^. 

Eq.dl) is similar to the generating functional obtained in Ref. ^ above threshold, except 
for new terms that we discuss in the next paragraph. The function C is the same as the 
mean field phase phase correlations as a function of time above threshold: C{vti — vt2) = 
{[(piti) — vti][(f){t2) —vt2]). Below threshold, the spatial average of the phase, {4>(t)), replaces 
vt. For any value of (0), the function Fmf{{4')) is the force that would result in an average 
displacement (0) in the mean field approximation to Eq.(|]), while the function F is the 
force that would produce the same average displacement in the full short-ranged model of 
Eq.(^. As mentioned in the previous paragraph, the difference between these two, which 
is the coefficient of the second term in Eq.(|^), can be found by the consistency condition 
($) = 0. This again is similar to the case above threshold. 

Despite the similarities, there are important differences between Eq.(|^) and the corre- 
sponding generating functional above threshold!. There is a non-zero 'mass term' in Eq.(^), 
r<l>$, which results in a finite response of to a spatially uniform dc force. Above thresh- 
old, where this term vanishes exactly due to time translational invariance of the system, 
xiQ = 0,ct;) has a l/iu singularity, corresponding to a finite response in the velocity of 0. 
The bare value of r is the inverse of the mean field polarizability x^(-F) (along the monoton- 
ically increasing path); r vanishes at threshold, although the manner in which it approaches 
zero at Ft depends on the distribution p{h) of the pinning strengths. The generating func- 
tional Z also differs from the form above threshold by various other terms that also vanish 
at threshold like <|)r$ does. For instance, the coefficient of the second term, F — Fmf{{4>)) 
has corrections of 0{[Ft — F^). These, however, do not affect the critical behaviour and are 
not shown in Eq.(^). Also, it is only when the mean field correlation and response functions 
are evaluated within the monotonically increasing approach to threshold that we obtain the 
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functional forms in Eq.(§) that resemble those above threshold. For other (non-monotonic) 
approaches, this will not be the case. Above threshold, where the steady state behaviour is 
unique, no particular path needs to be specified. 

The variable {4>{t)) in the arguments of the functions Fmf and C in Eq.(^) implicitly 
controls the dependence on F. It is necessary to increase adiabatically in order for 

the form in Eq.(||) to be valid; otherwise, there would be various transients that would have 
to be included in the equation. Since the increase of the force is much slower than all the 
dynamics of the CDW, the equal time response and correlation functions really give the static 
behaviour. Although it would seem that Eqs.(H) and (P) can be used to obtain responses 
to time dependent perturbations, there is an implicit constraint that the perturbing force 
must satisfy for our entire approach to be valid: e(ti) > ^(^2) for ti > t2- This constraint 
prevents us from obtaining the ac response, xi^y 

We now carry out a renormalization group analysis of Eq.(^ along the lines of Ref. |^. 
We rescale space and time as a; = bx' and t = b^t', and integrate out modes of all frequencies 
in the momentum shell between the (rescaled) momenta bA and A, where A is the upper 
cutoff in momentum. The rescaling of the fields $ and •I' are fixed in mean field theory 
by the requirement that all quadratic terms in the exponential have to be invariant under 
this transformation, except the <|>r$ term; the variation of this last term is interpreted as 
a change in F under renormalization. Fixing the other quadratic terms in the exponential 
yields $ = b^^^'^^'^^' and $ = 6^^'^/^$', with the dynamic exponent z fixed at 2 by comparing 
the scaling of the dt and terms. The consistency of the mean field scaling is verified by 
ascertaining that all higher order operators are irrelevant; this is indeed the case for d > 4, 
which is therefore the upper critical dimension. This upper critical dimension is the same 
as for F > F^. 

At d = A, the entire function C[(</)(ti)) + — [ti ^ t2)] becomes marginal. For d < 4, 
the critical behaviour is changed from its mean field form. As discussed in Ref. ^ for d < 4 
the function C splits up into a constant part Cs that controls the static distortions in the 
phase, and a functional operator C(0) that controls the dynamics. The two of these decouple 
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from each other, and scale differently. For the dynamics above threshold, where the steady 
state solution is periodic, it is necessary to preserve the invariance under the transformation 
(j) ^ (j) + 271. For d < dc, the upper critical dimension, fluctuations in (p scale in the same 
way as (0) = vt; this then requires $ to be invariant under rescaling, as discussed in Ref. 
^. The operator whose coefficient is the constant part of C, Cs, is then relevant, leading to 
anomalous scaling of the static correlations. 

Below threshold, it is precisely these static correlations that we are interested in. Since 
the mean phase (0) reaches a constant value for any fixed value of F < Ft in steady state 
(with {(p) increasing with F), the periodicity of the phase variable need not be preserved 
under the rescaling of the renormalization group. Accordingly, we choose a scaling for the 
fields $ and $ different from that of Ref. ^ 

I. = b-'^^^^'^', $ = 62-d/2$', (10) 

The scaling of the $ field is fixed by requiring the invariance of the Cg term, which is 
unrenormalized by loop corrections. (This is because all loop terms involve differences 
between the function C at two different values of its argument, which are not changed by 
an addition of a constant Cg- Although it is conceivable that singular corrections could 
affect this result, we assume here that this is not the case.) As discussed in Ref. ^ any 
renormalization of the term from the operator C (as well as from higher order irrelevant 
operators) must be of the form ^dt^, and thus the $V^$ term also has no loop corrections, 
fixing the scaling of $. The scaling of the term yields a renormalized polarizability 
x' = x/b"^- Finally, the {F — Fmf)^ term, which does receive loop corrections, is relevant, 
so that the scaling of F can also be obtained by power counting: [F — Ft)' = b'^^'^{F — Ft), 
thus yielding the correlation length exponent 

u = 2/d (11) 

and, from the scaling of x, 

7 = (12) 
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as mentioned in Eqs. and (0). 

Although it is possible to obtain these exponents directly from a very simplified model 
of CDWs0, the exponent z cannot be obtained from a similar analysis and, as we shall now 



see, is nontrivial. A more detailed discussion of Ref. |TT] is given in Section V. The exponent z 
is not affected by the difference between the scaling here and that of Ref. ^ although under 
the scaling used here, the operator C is dangerously irrelevant, flowing towards a singular 
function of zero amplitude, with constant second derivative C"{(f)). (With the scaling of 
Ref. ^, C flows towards a regular function of constant amplitude, with the same second 
derivative.) The 0(e) result for z, obtained in Ref. ||, depends on the second derivative of 
C, and is given by 

z = z{F > Ft) = 2-e/3 + 0{e'^). (13) 

The breakdown of two sided scaling for CDWs, with u = 2/d below threshold, and 
u = 1/2 above threshold, is now seen as the result of the different static and dynamic 
behaviour of the system, rather than being due to any fundamental difference between the 
properties above and below threshold. Even above threshold, the static correlation length 
still exists, and was indirectly obtained in Ref. ^ as the length that controls the flnite size 
scaling behaviour. The static correlation function above threshold varies as a simple power 
law with distance, and does not show any crossover at the static (or dynamic) correlation 
length, so that the static correlation length cannot be obtained directly. 



III. NUMERICAL RESULTS 

We now compare these results with those from numerical simulations. Previous simula- 
tions have often been carried out with the continuous variable dynamics of Eq.(|^), with x 
discretized to a lattice. However, it is possible to obtain much more accurate results with 
automaton models, which may be viewed as a singular limit of Eq.(^, in which the pinning 
potential at any site has very narrow and steep wells, so that (p{x;t) = P{x) + 27rm(a;;t), 
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with m an integer variable. A further approximation is made by discretizing time; B the 
dynamics of m{x] t) is then 

m{x; t) = m{x; t - 1) + 1 if F + V^[27rm(x; t - 1) + f3{x)] + h{x) > 

= m{x;t-l) ifF + V^[27rm(x;t- 1) + < 0. (14) 

Here the operator is the discrete Laplacian. This automaton model (and variants of it) 
are numerically much more efficient than direct simulations of Eq.(^. Since for small d, 
the critical behaviour is independent of the form of the pinning potential used, we shall rely 
mostly on the automaton model simulations in this section, referring to continuous dynamics 
only for comparison. 

It is more convenient to express the dynamics in terms of a local "curvature" variable 
c{x)= ^ [m(y) -m(x)] +Int[{F + V^/?(x) + /i(x)}/27r] (15) 

ye(xy) 

where the sum over y is restricted to sites neighbouring x. The first term in this equation 
is the discretized form of V^m. By construction, c{x) is an integer valued variable. From 
Eq. (p!4D , we see that the dynamics can be expressed in terms of c(x) as 

c(x) c{x) — n 

c{y) c{y) + iyy e {xy) if c{x) > 0. (16) 

Here n is the number of neighbours that the site x has. We can view this as a sandpile 
model0 if we treat c{x) as the height of sand at the site x; when this exceeds a preset 
threshold (zero in Eq. ([T6|) ) one grain of sand falls off from x to each of its neighbouring 
sites. As F is increased, Int[(-F + V^/3(x) + /i(x))/27r] increases by 1 whenever its argument 
crosses an integer. From Eq. (pTsP , we see that increasing F towards threshold is equivalent 
to dropping grains of sand randomly on different sites of the system. There is, however, a 
slight difference in the manner in which sand is added as compared to standard sandpile 
models: since the randomness from V^/?(x) + h{x) is quenched, sand grains are dropped in 
a random order on the different sites of the lattice one by one, with no site being revisited 
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twice in a cycle. When every site has been visited exactly once, the whole cycle is repeated 
in the same order. In contrast to this, in standard sandpile models the addition of each 
sand grain is an independent random process, and there are no correlations between the 
sites chosen. For a sufficiently large system, where the repeating cycle is very long, it is 
reasonable to expect the difference between cyclic and uncorrected addition of sand grains 
to be inconsequential .0 

Another difference between CDWs and usual sandpiles is that open boundary conditions 
are normally used for sandpiles, whereas periodic boundary conditions are used for CDWs. 
With periodic boundary conditions, a grain of sand that exits the system at one boundary 
reenters at the opposite boundary. Again, when the system size is much larger than the 
characteristic avalanche size, this difference should not affect the behaviour .BhI In this 
section, we shall use numerical results from CDW and sandpile models, both with periodic 
boundary conditions, to compare with our analytical results. Comparisons with earlier 
simulationJiiB on sandpiles with open boundary conditions will be discussed in subsequent 
sections. 

Figure 1 shows a finite size scaling plot for the polarization P{F, L) (measured from the 
initial state m = at F = 0) for c? = 2 as a function of the driving force F and the linear 
size of the system, L. The finite size scaling of the polarization should be of the form 

p = i/r^+ip(L|/r) (17) 

where / is the reduced force, F/ Ft — 1. The collapse to a scaling form in the numerical data 
shown in Figure 1 is not very good. This is probably due to the fact that the polarization in 
any state involves the polarizability integrated from the initial to the present state. Correc- 
tions to scaling might thus be expected to persist in the polarization until the system is very 
close to threshold. Figure 2 shows a similar scaling plot for the polarizability, obtained by 
a numerical derivative of the polarization. The data here scale better when L\fY is large, 
but show a large scatter very close to threshold due to numerical uncertainties. Figure 3 
shows a scaling plot for the polarization in the sandpile version of the dynamics of Eq.(|l^); 
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the scaling form works much better here. From the data in Figure 3, we obtain estimates 
for the critical exponents: ^{(1 = 2) = 1.98 ± 0.03 and ^[d = 2) = 0.98 ± 0.03, in agreement 



with Eqs.(|l2|) and (0) 



Simulations with the continuous dynamics of Eq. have larger numerical uncertainties. 
As detailed in Ref. ^, a log-log plot of the polarization as a function of / appears to yield an 
exponent 7 = 1.8 ± 0.15, although the slight upward curvature seen as F ^ Ft allows for 
the possibility that the asymptotic value of 7 is indeed 2, as given by Eq.([l2|). The collapse 
of the data to a single scaling function is also not very satisfactory, but can be used to obtain 
the correlation length exponent u as 1.0 ±0.1. A more accurate value for u is obtained from 
the width of the distribution of threshold fields F^ for various systems of size L with different 
realizations of randomness, since one expects finite size effects to become important when 
the characteristic size of an avalanche is of the order of the size of the system. Fitting AF^ 
to the form yields u = 1.01 ± 0.03.1 The results are thus consistent with Eq. (p!T|) . 

For three dimensional systems, the numerical uncertainties in the behaviour of the po- 
larizability are larger. Here we show only the results for simulations on the sandpile model 
(with periodic boundary conditions) for which, as in two dimensions, the scaling is better. 
Figure 4 shows a scaling plot of the polarizability, from which we obtain 7 = 1.31 ± 0.08 
and u = 0.68 ± 0.04. Since it is possible to determine the location of the threshold field very 
accurately, one can also plot the polarization at threshold as a function of system size. As 
shown in Figure 5, Ps{FxiL) appears to scale as L^-^"^, consistent with Eqs. (|ll|) and (p!2D 



for d = 3. (The error bars shown in Figure 4 are from the statistical scatter in the data; 
owing to the small range in L, this is too conservative an estimate for the true uncertainties.) 
The exponent v can be obtained independently from the scatter in the threshold field to be 
V = 0.68 ± 0.02, as shown in Figure 6. 

From the duration of an avalanche as a function of its size, we obtain the dynamic 
exponent z. The size s of an avalanche, defined as the total number of sites participating in it, 
scales as l'^, where / is its linear extent, while z is defined through t{l) ~ l^, so that t{s) ~ s^^'^. 
Figures 7(a) and 7(b) are log-log plots of t{s) versus s for CDWs in 2 and 3 dimensions 
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respectively. From these we obtain z{d = 2) = 1.32 ± 0.04 and z{d = 3) = 1.65 ± 0.06. This 
agrees very well with our result z{F < Ft) = z{F > Ft); z{F > Ft) can be found from the 
velocity exponent (3 through z{F > Ft) = 2(3, and (3 is known independently from numerical 
simulations above threshold to be (3{d = 2) = 0.64 ± 0.03 and (3{d = 3) = 0.81 ± 0.03.0^1 



IV. DISTRIBUTION OF AVALANCHES 

We now turn to the distribution of avalanches of different sizes in the critical region. As 
discussed earlier, there are differences in the forms for the continuous variable and automaton 
models. This is because the total increase in the phase at any site is bounded above by 27r 
in the continuous variable models. No such constraint exists for automata, for which several 
avalanches of roughly the same size are grouped together into a single large avalanche. This 
may lead to different number distributions for avalanches in the continuous and automaton 
models. We first deal with the avalanche distribution for automaton models, returning later 
in this section to continuous variable models. 

We conjecture a scaling form for the number distribution of avalanche^: 

1 r/s 

^(^5/M^ = ^^(^i/r')v ^^^^ 

where f)dsdf is the total number of avalanches of size between s and s + ds that occur 
in a unit volume of the system when the reduced force {F/Ft — 1) is changed from / to 
f + df. This equation defines the exponent n. (Here we have assumed that the same u is 
associated with the characteristic size of the avalanches as well as the finite size effects.) 

We define the moment of an avalanche as the total change in polarization that results from 
it. This is related to the change in the phase at every site in the avalanche by 6P = J2i=i{^'Pi)- 
For continuous dynamics, where 6(f)i < 2tt for all the sites, the moment scales as the size of 
the avalanche. For automaton models, the two are related by a nontrivial exponent: 

6P ~ s^. (19) 



18 



Equations (|1^) and (|T9|) can be used to obtain an exponent identity relating k and F to 
7. The total change in polarization upon changing the force by df is found by integrating 
Eq.(|18]) to be of the form AP ~ Comparing with Eq.(^, we obtain! 

^={Td-K)u. (20) 

Substituting Eqs. ([l^) and (|ll|) in this equation yields 

Td-K = 2. (21) 

Figure 8(a) shows a finite size scaling plot of the distribution of avalanches for the 
automaton model in two dimensions. By fitting to the form of Eq. (pISD , we obtain 

K(rf = 2) =0.36±0.03. (22) 

Figure 8(b) shows a similar plot for d = 3, from which we can obtain 

k(c^ = 3) = 1.00±0.06. (23) 

The values of u that are obtained from these scaling plots are i>{d = 2) = 0.98 ± 0.03 and 
z/((i = 3) = 0.68 ± 0.03, consistent with Eq.([TT]) (as well as the scatter in Ft). The scaling 
of the moment of an avalanche with its size can be used to obtain numerically 

r(rf = 2) = 1.15±0.05 (24) 

and 

r(rf = 3) = 1.00±0.02 (25) 

which, together with the numerical values for k, are consistent with Eq.(plD. In the limit as 
F —>■ Ft, a power law distribution for the avalanches is obtained. If n{x) is not a singular 
function of its argument x for small x, the limiting form of the distribution can be used to 
obtain k and F. If, however, n(x) scales as x^"~^^^'^'^'^ for small x, the apparent power law in 
the avalanche distribution would change from k to k — (« — l)/2z/, while leaving F unchanged. 
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A singular form for n(0) would imply that the total number of avalanches generated by a 
small increase in F (which is dominated by small avalanches) would be singular as Ft is 
approached; we have verified numerically that this is not the case for the automaton models. 

It is interesting to note that Eq.(^) has been derived earlier by Majumdar and DharEl for 
sandpile models, by a completely different method that relies only on properties of clusters at 
the critical point .il Since sandpiles at the critical point are current driven, J = 0+, instead 
of being tuned by a driving force F, any prefactor in n{s) that is singular in / would not 
be seen in the scaling. The fact that Eq.(^ll) and the result of Majumdar and Dhar are 
identical can thus be taken as a proof that n{0) is not singular. Our results in Eqs. (|22D 
through (p5D are also consistent with numerical simulations on sandpile models in two and 
three dimensions. 

As mentioned in the previous section, the dynamic exponent z can be found numerically 
from the scaling of the duration of an avalanche with its size; the results in both two and 
three dimensions are fairly close to a first order truncation of the e-expansion result of 
Eq. (pi^D . Even at the critical point, the duration of the avalanches can be used to obtain 
z; this allows us to compare with the result for two-dimensional sandpilesEl that z = 5/4, 
which agrees with our numerical results. For three dimensions, there is no exact calculation 
of z at present to compare our numerical results with. However, previous simulations on (a 
transformed version of) the sandpile model with open boundary conditions^ have yielded 
z{d = 3) = 1.62 ± 0.01, in agreement with our results. 

For continuous variable models, we have seen that the bound on the change in phase at 
any site in an avalanche constrains Fc to being 1. (The subscript c here denotes continuous 
variable dynamics.) From Eq. (pT]) , we then obtain Kc = d — 2, so that both the exponents 
involved in the avalanche distribution are determined.^ This result is only true if the driv- 
ing force is increased infinitesimally . In practice, the force is increased in small steps; no 
matter how small the step size, one presumably eventually crosses over into a regime where 
successive avalanches at the same site trigger within a single force step, and the exponents 
change to those of the automaton models. Although the analytical treatment of Section III 
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is applicable only for continuous variable dynamics, a proper analysis of the singular part 
of the non-linear response to finite increases in F might yield the automaton values for the 
exponent F; we have so far been unable to carry this out. Note that the numerical results 



fact, in three dimensions they are consistent to within the numerical uncertainties. 

Instead of triggering avalanches by increasing the driving force, it is possible to do so 
by applying small random 'kicks' to the phases at different sites, which may be considered 
to mimic the effect of thermal noise. Remarkably, for small thermal noise, the character- 
istic linear extent of an avalanche still scales with the distance from the zero temperature 
threshold force, but with a different correlation length exponent: 



with 7^ z^- The 'thermal' correlation length exponent has been found numerically by 



MyersQ to be ~ 0.50 in two dimensions. As shown in Figure 9, the same result is true 
for (i = 3 as well. A proper understanding of this phenomenon is still lacking; however, 
it is possible that considerations about the equilibrium behaviour (at finite temperature) 
in the absence of a driving force may be applicable here. At low temperatures, one might 
expect the presence of a driving force to produce an overall forward drift of the CDW, 
while still preserving the equilibrium properties locally. In equilibrium, the dynamics can 
be expressed in terms of correlations in the pinning potential. This is in contrast to our 
zero temperature analysis, where the physically relevant correlation function C is related 
to force-force correlations. The potential-potential correlation function is related to the 
function C integrated twice; as shown in Ref. in equilibrium it is necessary for this to 
be periodic in its argument. Changing C by a constant Cg is therefore no longer allowed in 
equilibrium. Setting Cg to zero implies that the static correlations no longer scale differently 
from the dynamic correlations, and = 1/2. Further work on this, especially on crossovers 
from equilibrium to nonequilibrium properties, is clearly required. 



in Eqs.(^) through ( p5D are not very far from the trivial values F = 1 and k = d — 2; in 



Cth ~ {Ft - F) 



(26) 
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V. DISCUSSION 



Most of the work on sandpile models has been on the properties at the critical point. 
However, Tang and Baklll have considered sandpile models away from the critical point; 
through a combination of scaling arguments and numerical simulations, they have obtained 
results similar to some of ours. For instance, by considering the diffusive nature by which 
sand grains propagate through the system, they obtain the exponent identity 7/z/ = 2, 
which is the same as from Eqs.(|T^ and (|TT|). With the assumption that the distribution 
of avalanches is independent of system size for sufficiently large systems, and the implicit 
assumption that the moment of an avalanche is proportional to its size, they obtain an 
exponent identity connecting 7 and k, which is a special case of the general result derived 
in the previous section. The first of these two assumptions is equivalent to h{0) not being 
singular; we have seen in the previous section that this is correct, both by comparing our 
scaling law to that of Majumdar and DharEl and by direct numerical verification. The second 
assumption, that F = 1, is not correct, even though it is approximately valid for d = 2 which 
is the case Tang and Bak focus on (and actually much better for d = 3). A few other results 
that they obtain are for the behaviour of transients, and cannot be compared with our work. 
The authors have also noted that if the critical behaviour in these nonequilibrium models 
were to be like that of equilibrium statistical mechanical systems, various exponent identities 
relating f3 to other quantities would result. These exponent identities can be compared with 
our results, and are found not to be valid. 

Surprisingly, the value obtained by them for u through numerical simulations is uld = 
2) = 0.7, which disagrees completely with our analytical and numerical result z/((i = 2) = 1. 
It is not clear what this difference is due to, although a similar value oi v {v = 0.76) has 
apparently been observed in other simulations as well in which the sandpile is perturbed 
in a manner different from ourS. As discussed in the previous section, with different 
types of kicks on the sandpile (that may or may not have physically meaningful analogs 
in the continuous variable models) v can be very different; this might be the source of the 
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discrepancy. Our simulations on sandpile models, detailed in the previous section, yield 
u{d = 2) ^ 1. 

A related problem to CDWs, that of the depinning of interfaces in random media, has 
been analyzed in the moving phase by methods similar to those used for CDWs.S'S Unlike 
the case for CDWs, the statics and dynamics do not behave differently; physically, this is 
because as the interface moves forward it experiences different regions of impurities and is 
thus unable to build up anomalous static correlations. Formally, this is seen by the absence 
of an unstable direction at the fixed point of C when the appropriate boundary conditions 
are taken. Since this separate scaling of static and dynamic correlations is the source of the 
difference between z/(F > Ft) and z/(-F < F^) for CDWs, two-sided scaling will be valid for 
interfaces. The relation 7/1/ = 2 can be verified to be still true by the procedure in this 
paper. 

Parisi and PietroneroS have considered a model for CDWs which is equivalent to 
choosing the pinning potentials V{(f)) at different sites to be of the form V{(j); x) = 
h{x)(j) — 2TTh{x) J2n ^(0 " 27rn), so that the phase at any site experiences a constant 

retarding force from the pinning potential, equal to —h{x). In addition to this, there is the 
elastic force from the neighbouring sites, and the external driving force. (Although they 
have a ^-function constraint that prevents the total force at any site from being negative, we 
have seen that when threshold is approached monotonically the phases at all the sites always 
move forward, so that this constraint is inconsequential for the critical behaviour.) Although 
this model correctly yields 7 = 4/(i and u = 2/d, the dynamic exponent z is wrong, being 
trivially equal to 2. The reason for this partially correct result can now be understood: after 
averaging over the randomness, their model corresponds to the one we have considered here, 
with the dynamic part of the correlation function, C{(f)), set to zero, but with Cs 7^ 0. This 
does constitute a fixed point, at which loop corrections from C{(f)) are absent, so that z = 2. 
However, this fixed point is unstable to the one that we have analyzed in this paper, so that 
the model that they have considered constitutes a pathological limit for pinning potentials. 
In particular, it is possible to replace the vertical drop in V{(f)) at = /5 + 27m with a 
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very steep section of non-zero width. The critical behaviour of this model has been found 
numerically! to the same as for the other models that we have considered here, although the 
width of the critical region above threshold can be seen through a perturbation expansion 
to vanish as the width of the steep sections grows smaller.0 

As mentioned in Section III, due to the restriction that the driving force must be mono- 
tonically increasing in time, it is not possible to obtain the ac response xi^^ ^) from our 
analysis. It is not clear at present if some variation of our approach may prove adequate 
for this purpose. It is important to remember that the universality that exists for various 
different models with regard to the monotonic approach to threshold is considerably limited 
for xi^j ^) even in low dimensions. 

In this paper, we have analyzed the critical behaviour of charge density waves as the 
depinning threshold is approached monotonically from below. Our success in this was due 
to the fact that, for this special path to threshold, the existence of many metastable states 
was argued not to invalidate the perturbation expansion we carry out. Although for any 
particular problem one must consider afresh whether such a perturbation expansion is valid, 
it may be possible to use similar techniques for other systems. 
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FIGURES 

FIG. 1. Attempted finite-size scaling of the polarization P for the CDW automaton in cZ = 2 

witli best fits 7 ~ 2 and u k, 1. No set of exponents gives a single curve, due to the large finite size 
effects and the arbitrary constant in the definition of the polarization. Symbols indicate the size 
of the system, while the numbers in parentheses indicate the number of samples averaged over. 

FIG. 2. Scaling of the polarizability (numerical derivative of the data in Fig. 1) in the 

d = 2 CDW automaton model with fitted exponents 7 = 2 zb 0.15 and f = 1 zt 0.1. The scatter at 
low scaled fields L\f\'^ is due to statistical fluctuations. 

FIG. 3. Scaling of the polarization in the d = 2 sandpile automaton model. The scaling 
collapse is much better than in the d = 2 CDW model. From the range of exponents for which the 
scaling collapse is within finite-size errors and statistical uncertainties, we find 7 = 1.98 ±0.03 and 
v = 0.98 =b 0.03, in agreement with the analytical results in the text. 

FIG. 4. Scaling of the polarizability in the d = 3 sandpile automaton model. From the range 
of exponents for which the scaling collapse is within finite-size errors and statistical uncertainties, 
we find 7 = 1.31 lb 0.08 and u = 0.68 ± 0.04, in agreement with the analytical results of the text. 

FIG. 5. Polarization at threshold for the d = 3 sandpile automaton with periodic boundary 
conditions as a function of linear size L. Error bars represent la statistical uncertainties. The fitted 
slope in the figure indicates only statistical errors and are an underestimate of the true uncertainty 
in the ratio (7 — l)/z^. From such data, we estimate (7 — l)/i/ = 0.57 it 0.08. 

FIG. 6. Relative sample-to-sample fluctuations AFt/Ft in the threshold held Ft as a function 
of linear size for the d = 3 sandpile model with periodic boundary conditions. From the fltted slope, 
the finite-size correlation length exponent 0.68 it 0.02 is obtained. 
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FIG. 7. Numerical derivatives of the logarithm of the duration t of an avalanche with respect 
to the logarithm of the volume s, for (a) the d = 2 CDW model and (b) the d = 3 sandpile model. 

The approach to a constant at large s gives the scaling t ~ s^/'^, with z = 1.32 it 0.04 in d = 2 and 
z = 1.65 ± 0.06 in d = 3. 

FIG. 8. Scaled avalanche distributions seen at different fields for (a) 100 samples of size 512^ 
in the d = 2 CDW model and (b) 110 samples of size 128^ in the d = 3 sandpile model. Avalanches 
result from the (adiabatically) increasing field. We find n = 0.36 ± 0.03 and u = 0.98 ± 0.03 for 
d = 2 and «; = 1.00 ± 0.06 and u = 0.68 ± 0.03 for d = 3. 

FIG. 9. Scaled distribution of avalanches for the d = 3 sandpile automaton model at fixed field 
(total height). Avalanches are in response to individual "thermal kicks" which increase the phase 
at one location (by redistributing the curvature variable). The scaling collapse is within errors for 
K = 0.99 ± 0.03 and i/ = 0.50 ± 0.02. 
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Fig. 7a - Narayan and Middleton 
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Fig. 7b - Narayan and Middleton 
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